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Abstract. A model Hamiltonian representing the Cu spins in La2Cu04 in its low- 
temperatm'e body-centred orthorhombic phase, that includes both spin-orbit generated 
Dzyaloshinskii-Moriya interactions and interplanar exchange, is examined within the 
RPA utilizing a Tyablikov decoupling of various high-order Green's functions. The 
magnetic susceptibility is evaluated as a function of temperature and the parameters 
quantifying these interactions, and compared to recently obtained experimental data 
of Lavrov, Ando and collaborators. An effective Hamiltonian corresponding to a 
simple tetragonal structure is shown to reproduce both the magnon spectra and the 
susceptibility of the more complicated body-centred orthorhombic model. 

PACS numbers: 75.25. -hz, 74.72.-h, 74.72.Dn, 75.30.Cr 
1. Introduction 

Experimental studies of the magnetic and electronic properties of the cuprates continue 
to produce new and unexpected results that spur on theorists in their attempts to 
understand and describe the underlying orderings and excitations that may be involved 
in the pairing leading to high-temperature superconductivity. One experiment, that 
was the motivation for the work that we present in this paper, concerned the (zero- 
field) magnetic susceptibility of undoped La2Cu04 - it was found [1] that the magnetic 
response of undoped La2Cu04 was highly anisotropic, and that this anisotropy persisted 
well above the Neel ordering temperature. Further, they found that this anisotropy 
persisted in the weakly doped state. 

The importance of this result can be recognized if one notes the ongoing efforts of 
various researchers in understanding the origin and nature of so-called stripe correlations 
that are found in some cuprates (for a recent review of this problem, see reference [2]). 
That is, if the undoped state has a highly anisotropic magnetic susceptibihty, can it 
really be that surprising that "spin stripes" are also present when the system is doped, 
and if not, what role does the anisotropic magnetic response play in the formation of 
stripe-like structures? 
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Previously, we have examined [3] the origin of this magnetic anisotropy by 
considering a single Cu02 plane utilizing a magnetic Hamiltonian that contains spin- 
orbit generated Dzyaloshinskii-Moriya (DM) interactions [4, 5], and near- neighbour 
superexchangc. If one includes both the symmetric and anti-symmetric DM interactions 
one finds that a true phase transition (at a non-zero temperature) occurs to an 
antiferromagnetic (AFM) state, wherein the AFM moment lies in the plane, with a weak 
parasitic ferromagnetic moment generated by a small canting of the moments out of the 
plane. Within mean-field theory, linear spin-wave theory, and within the RPA utilizing 
the Tyablikov decoupling scheme, we determined the magnetic susceptibility, and found 
(i) that it was indeed highly anisotropic, even when the DM interactions were small 
compared to near-neighbour intraplanar superexchangc, and (ii) quantum fluctuations 
produced a substantial modification of the susceptibility as one used a more and more 
"sophisticated" theory [3]. Other potentially important terms {e.g., cyclic ring exchange 
[6]) that could have been included in that paper are discussed at the end of this paper. 

In this report we focus on an augmented model that now includes the third 
dimension and the full body-centred orthorhombic structure of La2Cu04. Our 
motivations for doing so are as the following. I - Although one can produce a true phase 
transition within a model that accounts for only a single plane, the interplanar exchange 
interactions can also produce a phase transition in approximately the same temperature 
range. That is, some researchers have suggested that both the (unfrustratcd) interplanar 
exchange and the DM interactions arc of comparable strength, and thus there is no good 
reason to exclude cither of these terms in our model Hamiltonian (see [7] and references 
therein). II - As we discuss below, there are two different "near" -neighbour interplanar 
exchange constants, and these are different in different directions. Thus, this difference 
will be a source of magnetic anisotropy, and it is necessary to determine the extent of 
this anisotropy through a calculation that includes both the DM interactions and the 
interplanar exchange. Ill - As mentioned above, our previous work noted the strong 
effect of quantum fluctuations in a two-dimensional model. Since one expects such effects 
to be larger the lower the dimensionality of the system, it is possible that this behaviour 
is reduced in a full three-dimensional model. In this paper, again with the RPA utilizing 
the Tyablikov decoupling scheme, we have completed the requisite calculations for this 
more complicated but also more realistic magnetic Hamiltonian. 

Our paper is organized as follows. In the next section we summarize the formalism 
necessary to analyze this problem; although somewhat similar formalism is presented our 
previous paper [3], when going from 2D to 3D the analysis is much more complicated, 
and it is thus necessary to present the required equations that must be solved. (Some 
aspects of the calculations have been put into various appendices.) In the subsequent 
section we present the results of a detailed and exhaustive numerical study of the 
resulting formalism for reasonable parameter values. Then we suggest a simpler 
model Hamiltonian, one for a simple tetragonal structure which avoids the frustrated 
interplanar AFM interactions of the original body-centred orthorhombic structure. 
Finally we conclude the paper by discussing the key results that we have obtained, and 
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then provide a comparison between the predictions of our theory and the experiments 
of Lavrov, Ando and co-workers [1] . 

2. Model and Methods 

2.1. Model Hamiltonian 

We describe the magnetic structure of the La2Cu04 crystal in the low-temperature 
orthorhombic (LTO) phase by using an effective spin-| Hamiltonian for the Cu^+ 
magnetic ions of the Cu02 planes defined by 



In this equation Sj denotes a spin at site i, and sites labelled as ii and ji are in the "first" 
plane while ^2 and j2 are in the "second" (neighbouring) plane; the notation {ia,jp) 
refer to near-neighbour sites. This Hamiltonian is written within the xyz orthorhombic 
coordinate system shown in figure 1 (see right-hand side) and in figure 2(a), in what we 
refer to as the "initial representation" in the LTO phase. 

The various terms in the magnetic Hamiltonian given in equation (1) correspond to 
the following interactions. As was mentioned in the introduction, the orthorhombic 
distortion in the La2Cu04 crystal, together with the spin-orbit coupling, lead to 
the antisymmetric Dzyaloshinskii-Moriya (D term) and the symmetric pseudo-dipolar 
(r term) interactions within the each Cu02 plane [8, 9, 10]. These interactions together 
with superexchange one (J) can give rise to an ordered phase within a single Cu02 plane 
at some nonzero temperature [3] . In this long-range ordered state Cu spins are aligned 
antiferromagnetically in the ^/-direction, with a small canting out of the plane. Therefore, 
each Cu02 plane in a La2Cu04 crystal exhibits a net ferromagnetic moment, so- 
called weak ferromagnetism (WF) in the direction parallel to the c-axis of the Bmab 
space group (2;- axis in the initial coordinates). Due to the weak antiferromagnetic 
coupling between the planes, the net ferromagnetic moments of adjacent CUO2 planes 
are antiferromagnetically aligned and the system possesses no net moment. 

Each Cu spin has four sites above and below it in neighbouring planes. If all of 
these distances were equal the system would be frustrated because the ordering in one 
plane would not lift the degeneracy that would result in adjacent planes. However, in the 
LTO phase these distances are not all equal, and thus the interplanar coupling between 
nearest-neighbour spins depends on which pair of neighbouring sites are considered. 
That is, due to the small orthorhombic distortion (relative to the high-temperature 
body-centred tetragonal phase) some near neighbour sites arc closer together than other 
pairs (which are, technically, next-near-neighbour sites). In what follows we refer to 




(la) 



{lb) 
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the sites shown in figure 1, which allows for these ideas and the interplanar terms in 
equation (Ic) to be made clear. The distance between ji and j2 sites is slightly less than 
the distance between ji and 22, and thus the superexchange couplings are different, and 
in this paper we thus specify that neighbouring spins in the x — z plane {J± ) have a 
larger superexchange than do neighbouring spins in the y — z plane (J^): | J±| > | 
(see, for example, the discussion in reference [11]). As discussed in the introduction, 
and quantified in the next section, this difference immediately leads to an enhanced 
anisotropy of the magnetic susceptibility. 

We schematically illustrate the magnetic structure of the La2Cu04 crystal within 
the ordered state (T < Tn) in figure 1, where arrows represent the Cu spin structure; 
this ordered state is quadripartite, as we now explain. In our notations we label sites 
in a plane with the spin canting up as ii and ji, and correspondingly the sites of the 
nearest-neighbour planes with the spin canting down are labelled by indices ^2 and j2- 
In each plane sites with label i differ from the sites j by the spin orientation within 
the antiferromagnetic order. Clearly, the magnetic structure of the La2Cu04 crystal 
in the ground state can be represented by four different sublattices with different spin 
orientations, and in our calculations we will follow the notation that ii-sites belong to 
sublattice 1, ji-sites to sublattice 2, i2-sites to sublattice 3, and j2-sites to sublattice 4. 
The interaction of the spins from the sublattice 1 and 2 with the nearest-neighbour spins 
from sublattice 3 and 4, respectively, are described by term J±, and interaction with the 
spins from 4 and 3, respectively, are described by term J'j_. Each magnetic ion interacts 
with the four nearest neighbour sites within its plane and with eight ions (four above 
and four below) from neighbouring planes. 

To summarize, the above presented magnetic Hamiltonian describes the magnetic 
interactions within the each CUO2 plane in its first and second parts (equations (la, 16)), 
while the third part (equation (Ic)) takes into account the weak interplanar 
superexchange couplings. 

The structure of the Dzyaloshinskii-Moriya (DM) and the pseudo-dipolar 
interactions for the LTO phase are given by 



and 
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within the initial coordinate system [3]. The DM vector given in equation (2) alternates 
in sign on successive bonds in the a — b and in the a — c direction of each plane, as 
is represented schematically by the double arrows in figure 2(b). Thin arrows in this 
figure describe the in-plane antiferromagnetic order of the Cu spins, and are canted 
up/down from the in- plane order by a small angle. In the classical ground state of the 
LTO phase the absolute value of the canting angles are equal on all sites and are given 
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Figure 1. (Colour online) Magnetic structure of La2Cu04 crystal. Sites having 
different spin orientations are labelled by indices ii, ji in the plane with the WF 
moment in the positive z direction, and 12, 32 in the plane with the WF moment in 
the negative z direction. For each set of sites the a spin coordinate system within the 
characteristic representation (CR) is shown. The thin net is shown only to simplify 
the visualization of the canting in the spin structure. 



by the expression 

^ = itan-| , \. (4) 

2 \j + i(ri + r3)-|Ji/ ^' 

Following the scheme described in our earlier work [3], we perform rotations of 
the spin coordinate system in such a way that the new quantization axis (cr^) is in the 
direction of a classical moment characterizing the ground state. Hereinafter we will call 
such a representation as the "characteristic representation" (CR). Since four different 
types of spin orientations are present in the magnetic structure of the La2Cu04 crystal, 
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Figure 2. (a) Coordinates in the initial representation, (b) Thin arrows — the Cu 
spins, and open arrows — the DM vectors. 



we introduce four different spin coordinate systems (see the left-hand side of figure 1) 
given by the transformations equations (A.1-A.4) in Appendix A. Thus, each sublattice 
consists its own spin coordinate system. The model Hamiltonian in terms of these spin 
operators a in the CR is given by 



+ E |i(>^±+>^p)K^72 + ^n<2) + '\iJ'l--Jv){<^t - <^n'^l) + Jv<^n\ 

{h,j2) ^ ^ 



(«1.«2> {31,32) 

where we have used the following definitions: 

Jl — J3 rj Ja . .J1 + J3 /^v 

^ = 5 = (6) 

Jl = J + Ti, Jp = Ji cos 2^, (7) 

J,= £iZll+ (^j+£i±£i^cos2^ + Asin2^, (8) 

J3 = -£iZll+ (^j+£i±£i^cos2^ + -^sin2^, (9) 

J4 = -r2sin^+ ^cos^. (10) 
v2 
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The subscripts {i,j)ab and {i,j)ac in the summations of equation (5) imply the nearest 
neighbours in the ab and ac directions, respectively, as shown in figure 2(b). 

2.2. Mean field analysis 

In this subsection, wc present the results of the mean field approximation (MFA) for 
the above Hamiltonian by following the standard decoupling. That is, in the MFA in 
equation (1) we use 

- «) ^} + < (^?) - «) (^?), (11) 
where a and b can be equal to any of x, y, z. Then, the equation for the order parameter, 
to be denoted by r], within the MFA reads as 

rj = (a^) = Itanh |^Z[J2 + - J,]{a^)^ , (12) 

where Z — A is the in-plane coordination number, and (3 = 1/T. From this equation 
the Neel temperature at which rj vanishes can be written immediately as 

By applying a magnetic field sequentially in the a;, 7/, and z directions of each coordinate 
systems within the CR we can find the transverse and longitudinal components of 
susceptibility within all four sublattices. Using the relation between the components 
of susceptibility in the initial and characteristic representations given in equations (A.6- 
A.8), we obtain the final result for the zero-field uniform susceptibility within the MFA 
below the ordering temperature (T < T^^^) 

X''^^^=- r, (14) 

^ 4 Ji+J2+2Jx+(Ji-Jp)' ^ ^ 

y MFA ^ 1 sin^(^) ^ cos^^ sechmZri-J^fa} 

^ 4 J2- J3+2J±-2Jp 4 T + [J2 + J± + Jp] Sech^ {IZTjJrnfa}' 

zMFA_^ cos^W I sin'^ sech^{fZr/J^j4 

4J2+J3+2J^ 4 T- [J2-J±+Jp]sech2{fZr;J^y4' 
where we define 

Jmfa = J2 + 'J± — Jp , (17) 

and the equation for the order parameter rj given by equation (12). Wc used the "mfa" 
subscript in equation (17) because this combination determines the effective interaction, 
and thus the Neel temperature within the mean field theory (see equation (12)). 

The final results for the components of the susceptibility in the initial representation 
for high temperatures, that is above the ordering temperature {T > T^^-^), are 

X MFA _\ \ /-|o\ 

^ 4T + Ji + J^ + Ji' 

y MFA ^ 1 Sin^(^) , 1 COS^(^) Qq) 

^ AT-Js + J^-Jp AT+J2 + J± + Jp ^ ^ 

. MFA _ 1 cos^(g) 1 sin^(g) 

^ 4r+J3 + J± + Jp AT - J2 + J^- Jp ^ ' 
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In the limit T — > T^^-^ we obtain that the x component of the susceptibihty is 
continuous at the transition and is given by equation (14). The y component of the 
susceptibihty at the transition temperature reads 

MFA ^ 1 sin^(^) , cos^ e 1 

^ T^TMFA 4J2_ j3 + 2J^-2Jp 8 J^ + Ja' ^ ' 

Note that with respect to the pure 2D case [3] the z component of the susceptibihty 
does not diverge at the Neel point and also is continuous at the transition 

. MFA ^ 1 cos^(^) , 1 sin'(^) 

2.3. Random phase approximation 

In this part of the paper we use the technique of the double-time temperature dependent 
Green's functions within the framework of the random- phase approximation (RPA). In 
the imaginary-time formalism, the temperature dependent Green's function and the 
corresponding equation of motion for two Bose-type operators reads 

Gab{t) = {T^A{T)Bm, = 5{t){[AB]) + (T.[i7(r), A(r)]i?(0)), (23) 

where A(t) — e^'^Ae~^'^ is the operator in the Heisenberg representation for imaginary 
time argument r, and Tr is the time-ordering operator. 

By using the method proposed originally by Liu [12], we employ the perturbed 
Hamiltonian 

H( = HcR -fYl ' ^'i ^ sublattice 1 (24) 

to find the longitudinal components of the susceptibility in the CR. In the equation (24) 
/ is a small fictitious field applied to the spins of sublattice 1 only. In this paper we 
are studying the zero-field uniform magnetic susceptibility, therefore we restrict / to be 
constant and static. 

The Green's functions to be used in present calculations are 

G^ir) = (T.<(r)a7(0))^, G(r^(r) = (T.ar(r)a^7(0))/, 

G^,^.^(r) = (T.a+(r)a^7(0))/, Gj^^(r) = (T.aT, (r)a^7 (0))/, (25) 
= (T.<(r)a7(0))/, Gf-^(r) = (T.ar(r)a7(0))/, 

where {...)^ means that all expectation values are taken with respect to the perturbed 
Hamiltonian in equation (24). After an expansion in a power series of / the Green's 
function, e.g. G{^j^{t), reads 

GLi^) = Gf;]^{T) + (r) + Oif). (26) 



(0) 

GUi^) = G^,,,{r) + /GS;].^(r) + 0{f). (27) 



Since J-^ (r) = Gi^^j^ (r) , from now we drop the superscript and use 



dr 



Magnetic susceptibility of La2CuOi 9 
Also, we introduce 

K(r)K = «) + /Vn+0(/^), (28) 

where, due to the translation periodicity {al^) = rj, the order parameter at / = 0. 

Now let us find the equation of motion for the Green's function G{_^j_^{T). The 
equations for other functions can be found in the same way. Starting from the 
equation (23) we can write 

^^y^ = 26{T)6,,j,{aty + {TAHcKir), <(r)]a7(0))/ - (29) 

In order to solve this equation of motion we are following the RPA scheme, and using 
the so-called Tyabhkov's decoupling [13] which is given by 

{Tr<7!{T)at{T)ar{0)y - (<7r(r))^(r.<(r)<7-(0))^ = (<^f(r))^G'(,,(r). (30) 
After this decouphng is introduced, equation (29) is found to be 

" ^^Sm.jA^y-fGi^^ir) (31) 
E {2«(^))'Kn+5),. W - BGlr^+S)n(^)] + J2«+,(r))^G(,,(r)} 

Sab 

-E {2« + ^2«+,(r))^G(,,(r)} 

Sac 

-E {2«W)'[^G;,^.M-i^G£.(r)] - J,(.^(r))'GL,(r)} 

-E {-^M<ir)yGi,jT) + Jx(a4(r))/G(,.^(r)} , 

where Yls b ^^^^^^ ^o a summation over the nearest neighbours of the sites ii in the ab 
direction of the same Cu02 plane, and similarly for '^g^^ — see figure 2(b). Thus, in 
equation (31) all sites ii + 5 belong to the sublattice 2. The notation means sum 

over all sites ii^ from sublattice 3 which are nearest neighbours of sites ii that belong to 
the sublattice 1, and similarly for ■ 

Next, we perform the transformation into the momentum-frequency representation 
for the Green's functions and the spin operators: 

GUir) =^E^^2(fe,^^)e^'=-^^--^^-^^e-^-^ (32) 

k,m 

k,m k 

where the expansion of equation (28) and the linear response to the uniform perturbation 
expressed by vi(A;) = S{k)vi were taken into account. In the transformation given by 
equations (32, 33), the sum over k runs over ^A^ points of the first Brillouin zone, and 
LUn = 27rn//3 for n e Z are the Bose Matsubara frequencies. The equation of motion for 
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the Green's function G(^j^{t) in the momentum-frequency representation reads 

-icjmG(^{k, Um) = - 2ZAk[n + /vi]G'4(fc, Urn) + '^ZB^iv + fyiM^ik, Um) (34) 

- Qkiv + /vi]G'{2(fe, Um) + 2Z ibkiv + /vi]G'{2~(fc, o;^) 

- Z^Mf] + /V2] + J^lv + /V3] - Mv + fY,]]Gi,{k,u;J, 
where, as before, Z is the in-plane coordination number, and we introduce 









(35) 


ak = 


T' 4- T T' — T 


2 ^fc' 


(36) 


7fe = 


-(cos/c-r -|- COS /Cj/), ^fc = COS COS 


I 2 /' 


(37) 


l'k = 


- (cos kx — COS ky)^ Cfc = COS kz cos 


/ kx — ky\ 
\ 2 J' 


(38) 


write 


down the final equations for the zero-order in / Green's 


function 



Gi2{k,ujjn}, and the first-order one G^i2 i^i^m) 

icunt ^(1) ^ G12 f iiYl+'I±Yl-ipI±\G ^1 f _ f-Ii,-I±_:lp\\Q 

2Zr] ^'^ ~ 2Zr] \ 2 r] 2 r] 2 r] ) ^ r/ \2Z7] 1 2 2 2 J / 



.(1)- 
'■32 



+ {^"^^~^}^i2^ ^fc<^22^ ~ BkG 22 + akG42 — ibkG^42 ~ ^dkG^ 

(40) 

where in these equations we drop the wave vector and frequency dependencies for the 
Green's functions; that is G = G{k,Um) and G*-^-* = G^^\k,ujm)- 

In order to obtain a closed set of the equations for the zero and first order Green's 
function we should use the above described scheme for the all other functions in 
equation (25), and the final system of equations for the zero and first-order Green's 
function are given in the Appendix B in equations (B.1-B.3). The structure of the 
system for the zero-order functions is identical with the system of equations for the 
first-order ones, except for the free terms. Hence, the poles of the zero-order Green's 
functions (that determine the spectrum of the spin- wave excitations) G{k, cUm) are equal 
to the poles for the first-order ones G^^\k,u;jn), and are found to be 



£i,fe = 2Zr]uJi^k = \l tti,fc + , £2,fc = 2Zr)uj2,k = \/ai,fc - , (41) 

£3,fe = 2Zr]uJs,k = \J a2,k + \fP2^ , £4,k = 2Zr]L0i^k = \/ oi2,k - \/^%k , 

ai,k = o-l + {Ak-Jmfa/2f - {hk-duf - \Bu\\ (42) 

Oi2,k = 4 + (^fe+^m/a/2)' - {bk+dkf - |Sfe|^ 

= 4:[ak{Ak-Jmfal2) - {bk-dk)'^BkY - {2^Bkf[al - (6^-4)'], (43) 
/32,fe = 4:[ak{Ak+Jmfa/2) - {bk+dk)^Bkf - {2mkf[al - (6^+4)'], 
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within the notation of equations (35, 36), and the MFA-inspired definition of Jmfa — 
J2 + J± — Jp- 

The free terms in the first-order systems (see equation (B.3)) are determined by the 
zero-order Green's functions, and thus the first-order quantities G^^'' can be written down 
in terms of the solution for the zero-order system Appendix C, and the as-yet-unknown 
quantities Vi, V2, V3, and V4. 

To calculate Vi^2,3,4 we use a relation that connects v and the Green's functions 
G'W(fe,T = 0-), viz. ' ' 

-yi-j^T.^i\k,0-), Z= 1,2,3,4. (44) 

k 

After the substitution of the solutions of the systems of equations for the first-order 
Green's functions G^^\k,ujm) in Appendix D into the system of equations for v; in 
equation (44), the results are found to be 

Vi-V2-V3+V4= — 1 .„^n- r:.'^^,,/.. 0-^/7 , T^^r„r. (^5) 



where 



N ^ 

k 

4 



N 

k 







^Ee-^'^-°"A^n^m) 

m 


2rj- 


m 


m 


2rj- 


m 






— \Gi2 


2 1 \n- |2 \n- |2 \n |2 1 \ri |2 12 , \ri- 1 

+ 1*^221 "1*^121 ~|<^42| i-|Lf32| " |<J"42| +1*^321 


K^22^" 


-\Gi2\ 


2 \n- |2 1 \n- |2 \n |2 i i/^ 12, \ri- |2 \n- |2 

~l*-^22l +l*-^12l i-|<J-32| +|<J"42l ~l*-^32l 



(46) 



(47) 



k 



N ^ 

k 



\G22 


^+1^12!^+ ^22^+^^121^" 


^42^ — 


1/^ |2 1/^— |2 1/^— 1 
|^32| —\^A2\ ~l^32l 


\G22 


2 1 I/-Y |2 \ri- |2 \n- |2 
+ |Lri2| — |<J"22l ~l*^12l ~ 


^42^ — 


G^32^+ G42r+K^32^ 



(48) 

and all zero-order Green's functions G(fc, uim) are given in Appendix C. 

Now let us find the quantities which determine the linear response to a magnetic 
field applied to the one of sublattices {e.g., see equation (A. 5)). The longitudinal z 
components of the susceptibility in the characteristic representation are given by 



Xii 



df f=o df 



a" (J 



Xl3 = 



df /=o df 
where the expansion of equation (28) was used. The transverse x and y components of 
the susceptibility tensor are determined in the terms of Green's functions to be given 



= V2, (49) 

/=o 



= V4, 

/=0 
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by 

^ii'^"'^ ^ E /V<(r)<(0))dr, xrr'^"'= ^ E /Vr<(r)<(0))dr, (50) 

^la'^"'^ I E /V<(r)<(0))dr, Xi?"'- ^ E /Vr<(r)<(0))dr, 

where a = x,y. By substituting the solutions in Appendix C into the definition in 
equation (50) for the transverse components of susceptibihty, we obtain the result given 
in Appendix E. This result for the transverse components in the CR is exactly the same 
as the MFA calculations for the transverse components. 

Then, using equations (A.6)-(A.8) the components of the susceptibility in the initial 
coordinate system of equation (1) below the transition temperature are found to be 

= ^ (51) 

^ 4Ji+J2+2Ji+(Ji-Jp)' ^ ' 

^' = L.- J^'fa J, + 

These expressions for the components of susceptibihty include the as-yet-unknown 
value of the order parameter 77. It can be found directly; from the definition of the 
Green's functions we have 

Gnnir = 0-) = {a-al) = \-r], G^nir = 0") = ^ E'^22(fe,T = 0"). (54) 

k 

Substituting 6^22(^5 <^) from equation (C.2), and performing the summation on the 
Matsubara frequencies, the equation for the order parameter turns out to be 
1 14 / xi^fe \2n(£i,fe)+l / xi^k \2n(£2,fc)+l 

( X2,k \ 2n(e3,fc)-M / X2,k \ 2n(e4_fc) + l 1 
+ [y2,k+—?==] +[y2,k — 7=) \ , 

where 

xi,k = - 2afc[afe(Afe-J„/„/2) - {bk-dk)^Bk], yi,k = -(^fe-^m/a/2), 

X2,k = 2ak[ak{Ak+Jmfa/'^) - {bk+dk)^Bk\, y2,fe = (^fe+^'m/o/S), 

n(ei,k) = [exp(/3£i,fe) - 1^1,2, 3, 4. 

Since the order parameter {viz., the sublattice magnetization) is temperature dependent, 
it follows that the spectrum of elementary excitations (equation (41)) is also temperature 
dependent. 

The Neel temperature at which rj vanishes within the adopted RPA approximation 
is determined by 
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X2,k \ 1 ^ / X2,k \ 1 



By putting 77 — we find that the z-component of the susceptibihty 
equation (53) does not diverges at the Nccl temperature, whereas it diverges for the 
pure 2D model {J± = J'^ = 0). At the Ncel temperature all components of the 
susceptibility within the RPA are equal to the MFA results, the latter of which are 
given in equations (14, 21, 22). 

For completeness, we mention that the investigation of the model equation (1) 
within linear spin- wave (LSW) theory leads to the same structure of the susceptibility 
expressions as we found within the RPA in equations (51-53). The main difference 
between the results in RPA and LSW theory comes from the calculation of the 
longitudinal components of the susceptibility in the CR. The spin-wave theory gives 
unity in the denominator of the expressions in equations (45, 46), and S* = 1/2 instead 
of the order parameter 1] everywhere in the numerator Af^ and A/"^. The transverse 
components of the susceptibility in the CR are equal within the all of the MFA, RPA, 
and SW theories. 

When the temperature of the system is above the Neel temperature, Tjv, there 
still exists short-range antiferromagnetic order. To model such an order we follow 
reference [12] and introduce a fictitious field h pointing in the direction of the sublattice 
magnetization, that is the z direction in the characteristic representation. To this end, 
the Hamiltonian 

Hn^HcK-hY^at-hY^a^^ (57) 

i 3 

is used, and the limit /i ^ is taken after the calculation is carried out. 

Above the Neel temperature, we define a (different) order parameter by 

y ^\mi{2Zr]/h). (58) 
/i— >o 

By a procedure similar to the above presented [3] (that is, the RPA scheme below Tjv) 
we have found the equation for the order parameter and all components of the magnetic 
susceptibility in the paramagnetic phase. It is then possible to show that paramagnetic 
version of the equation for the order parameter in equation (55) leads to 

^-h^n (5.) 



2,fe 



y/K^^ ^Ik ^ 

X2,k \ 1 / X2,k \ 1 ] 

where in all expressions for xi,2, 2/1,2, P1.2 and cji_4 which determine equation (59) in 
the paramagnetic phase, we use a new definition for Jmfa (which we now call Jmfa) that 
reads 

Jmfa ^ J2 + J± - Jp + ^- (60) 
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The quantity y approaches infinity as the temperature is lowered to T^. Indeed, 
putting 3^ — > CX3 in equation (59) we find the temperature at which y diverges, which is 
identically equal to the Neel temperature. 

We have found (see below for numerical results) that for model equation (1) all 
components of susceptibility are continuous at the Neel point within the RPA. 

3. Numerical Results 

In this section we present the result of numerical calculations of the system modelled 
by the Hamiltonian of equation (1) based on the above- presented analytical formulae. 

3.1. Parameters regimes 

Firstly, let us consider the set of model parameters that appears in the Hamiltonian of 
equation (1), viz. in-plane parameters J, d and F, and out-of-plane parameters Jj^ and 
J'^. The in-plane parameter d that describes the antisymmetric DM interaction and 
parameters Fi 2,3 that give the pseudo-dipolar anisotropy, are of order 10^^ and 10^"^ 
respectively in units of J [9, 10], and it has been shown that the only combination from 
the pseudo-dipolar terms that affects the behaviour of the system is AF = Fi — F3 
[3, 14]. Thus, the in-plane part of the model, that is equations (la,16), can be 
completely described by the AFM Heisenberg model with the DM antisymmetric 
exchange interaction D and XY-like pseudo-dipolar anisotropy given by AF. 

In order to examine the behaviour of the system with respect to the out-of-plane 
parameters, we introduce the combination 

AJ^=J^-j; , (61) 

that describes the interplanar anisotropy interaction between nearest-neighbour spins 
which we refer to as the net interplanar coupling, and the combination 

Jx=J^+j; . (62) 

In our calculations we take A to be of the order 10~^ — 10~^ in units of J (see [7] and 
references therein). 

In this subsection we focus on the behaviour of order parameter 77, Neel temperature 
Tn, and susceptibility x with respect to the parameter within the RPA method 
(we present a detailed consideration of the dependence on AJ^ in a subsequent 
subsection). Firstly, we find that the order parameter and the Neel temperature are 
almost independent of the within a wide range of the model parameters. In figure 3 we 
show two representative plots for the order parameter and the susceptibility for certain 
values of the in-plane parameters. In each line of the figure 3a (that is solid, dotted, and 
dashed lines) we have simultaneously plotted five data sets, each with different values of 
the parameter that has been varied from zero up to 0.5 J ( J = 0, 0.01, 0.1, 0.2, 0.5). 
As one can see, for such a wide range of the parameter there is virtually no difference 
of the absolute values of the Neel temperature and the order parameter, whereas the 
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relatively small changes of the net interplanar coupling A in figure 3a strongly affects 
these quantities. 
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Figure 3. (Colour online) The (a) order parameter vs. T/J for the different values of 
AJ^: AJ_L=0- black solid line, AJ±=0.1 x 10" V - blue dotted line, AJ_l=0.4x10-3j 
- red dashed line, and each line consists of five data sets with Jj^ varying from zero 
up to 0.5 J. The (b) susceptibility, in units of 1/J, vs. T/Tm for AJ^=0, and for the 
following values of J^: J^=0 up to 0.02 — upper plot (many curves are superimposed 
on top of one another), J^=0.2 — middle plot, and J_l=0.5 — lower plot. In both 
figures (a) and (b) we have fixed d/J = 0.02 and AT/ J = 0.42 x 10"^ 

In case of the susceptibility, its dependence on Jx_ differs from that discussed above 
for the order parameter rj. Now, as is seen in figure 3b, the parameter J_|_ generates the 
constant shift in the ^"^^ components of susceptibility as well as the constant shift 
in the near the Neel temperature and within the paramagnetic region. However, for 
the reasonable values of the out-of-plane model parameters (that is J^, J'^< 10~^ J) the 
parameter J x. does not affect on the behaviour of the susceptibility with temperature. 

This latter result for the susceptibility can be shown in various limits from the 
above analytical results by taking into account the small magnitude of AJ^ and the 
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in-plane parameters d and AF with respect to J. In the zero temperature hmit one can 
write down the susceptibihty in the following form: 

11 11 d'^ 

^^'^°^4^7TX' ^^-°'^32{j-j^/4pAr + 2AJ/ ^^^^ 
Also, near the Neel temperature one finds 

At^Tn ~ 32 I J ^ j^i^y + 2AJx 4 2J + Ji' ^ ^ 

~i L_ + 1_L. (65) 

Almost everywhere within the above presented equations (63-65) one can ignore the 
contribution of J± < 10^^ J with respect to J; only the ^-component of the susceptibility 
at the Neel temperature is strongly affected by J±, as shown in equation (65). In 
fact, the upper plot in figure 3b consists of data sets with the different values of the 
parameter J± over the range of zero up to 0.02J, but these differing plots can not be 
distinguished from one another. 

Similarly, the parameter J± can be ignored in the expressions for the spin-wave 
gaps, as can be clearly seen from the following approximate formulae 

Ji 1 / . , _c[_ . , J| 



UJ 



l,fc^O 



nv2 



4fc^o ^\J-'^ K^r/2 + AJ^), (68) 



2 
J 



{j+:^}Ar/2. (69) 



Therefore, one can conclude that the affect of the parameter J± on the physics of the 
model is negligibly small. Consequently, without further trepidations the system can 
be studied using a fixed and representative value of this parameter, e.g. J± — 10"^ J, 
without having to be concerned with it changing our results. 

At the end of this subsection we discuss briefly the case of isotropic interplanar 
coupling {AJj_ = 0). In such a case the only difference in the susceptibility, with 
respect to a pure 2D model [3], is the finite value of at the Neel temperature. Then, 
only the in-plane anisotropics are responsible for the anisotropic magnetic properties 
of such a system, viz. the behaviour of susceptibility, order parameter, and spin-wave 
excitations. We emphasize the perhaps expected result, that for the 3D case with 
isotropic interplanar coupling, due to the frustration of interplanar coupling within the 
body-centred lattice, the effects of 2D quantum fluctuations and short-range correlations 
are very important, whereas the interplanar coupling is not. 
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3.2. Neel temperature and spin-wave excitations 

Now we present the results of our numerical investigations of the Neel temperature 
and the spin wave excitations and their dependence on the parameters of the in-plane 
anisotropics d/ J, AT/ J, and the out-of plane anisotropy AJj_/J. 

Figure 4a shows the Neel temperature obtained within the RPA scheme as a 
function of both AT/ J and AJj_/J. We found that the transition temperature Tjv 
depends on both the in-plane XY-like pseudo- dipolar anisotropy parameter AF and the 
out-of-plane anisotropy parameter AJj_, and changes of the same order (~ in 
AJ± and/or in AF produces considerable changes in the Neel temperature, T^. Further, 
the dependence of the Neel temperature on the DM parameter is shown in figure 4b: 
Tjv decreases as d increases for small d, but for larger values of the DM interaction, i.e. 
d » AJj_, AF, the Neel temperature Tjv increases nearly hnearly with d. The transition 
temperature into the long-range ordered state, TV, increases as the parameter AJ± 
increases since the net interplanar couphng that each spin feels favours to the AFM 
state. 

Figure 5 shows the zero-tcmpcraturc energy gaps in the long wavelength limit as 
a function of in- and out- of plane anisotropy parameters, and the resulting behaviours 
can be understood immediately from equations (66)-(69). Two modes, Si and £2, are 
almost independent of AF (see figure 5a, and equations (66,67)), but they show a strong 
dependency on the DM parameter, d, as seen in figure 5b. Since the canted angle goes 
as 9 (d/-\/2)/(2J), the modes £1, £2 are nearly linear in d. In the limit of zero DM 
interaction the mode £2 goes to zero and a Goldstone mode appears in the spin-wave 
spectrum, while the mode ei goes to a finite value, which is about 2Zri\/ JAJ^ (see 
equations (66,67)). Two other modes in the spectrum, £3, £4, are almost independent of 
the DM parameter of anisotropy, while they vary strongly with AF. In the limit AF = 
the mode £4 goes to zero and another Goldstone mode appears in the spectrum, while 
the mode £3 goes to the finite value of about IZrw/ JAV/2 (equation (68,69)). Since 
in the case of 3D model thermal fluctuations do not destroy the long-range ordering 
for T 7^ 0, the Neel temperature does not go to zero when a Goldstone mode £2 or £4 
appears in the spectrum. 

Lastly, we note that the plots of figure 5c show that two modes (£2 and £4) 
are independent of the net interplanar coupling AJ^, and two modes (£1 and £3) 
demonstrate a square-root dependency on AJ^. When the net interplanar coupling 
goes to zero, A J_l = 0, the two modes £1 and £2 become equal and describe the in-plane 
mode in the spin-wave excitations [15]. Similarly, the mode £3 coincides with £4 at 
AJx = and they correspond to the out-of-plane magnon mode [15]. 

3.3. Susceptibility 

Now we consider the temperature behaviour of the susceptibility and examine its 
dependency on different values of the in-plane and out-of plane anisotropy parameters. 
Our results for the y, and z components of the susceptibility within the different 
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Figure 4. (Colour online) The Neel temperature T/y, in units of J, vs. (a) AJ±/J 
- black solid line (fixed AT/ J = 0.42 x lO"^) and AT/ J - blue dotted line (fixed 
AJ±/J = 0.42 X 10"^), both for fixed d/J = 0.02, as well as vs. (b) the DM parameter 
d/J (fixed AT/ J = 0.42 x lO'^ and AJj_/J = 0.42 x IQ-^). 

approximation schemes, viz. RPA, LSW theory, and MFA, are presented in figure 6 
(T < Tat). We do not present the similar comparing for the x component of susceptibihty 
because of the pure transverse component (see equation (51)) has the same value 
within all mentioned approximations (below the transition temperature). On the 
other hand, the longitudinal (in the characteristic representation equations (A.7,A.8)) 
components of the susceptibility are involved in the equations for the components 
and x^ that leads to their different temperature behaviours within the different 
approximation methods (see below). 

Our results for the y component of the susceptibility, are shown in figure 6a. 
We find that at low temperatures the RPA analytical scheme, as was also found in pure 
2D case [3], is in good agreement with the linear spin- wave theory. Plots in figure 6a 
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Figure 5. (Colour online) The energy gaps, in units of J, as function of 

(a) XY-like anisotropy parameter AT / J {d/J = 0.02, AJ^/J = 0.42 x 10"^), 

(b) DM parameter d/J {AT/ J = 0.42 x 10"^^ AJ±/J = 0.42 x 10"^)^ ^nd 

(c) out-of-plane anisotropy parameter AJ±/J {d/J = 0.02, AT/ J = 0.42 x 10^^). 
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also show that RPA results agree with the MFA as one nears the transition temperature 
Tjv, and both RPA and MFA lead to the same magnitude of susceptibility at the Neel 
temperature. (It is worth noting here that transition temperature T^r within the MFA 
approach, where Tj^ = J2 + J± — Jp ^ J, is almost independent of the anisotropy, in 
contrast to the RPA scheme where T/v is very sensitive to the anisotropy parameters (see 
figure 4).) One can also see that in the zero temperature limit all approximations used 
in the paper go to the same value of susceptibility approximately given by equation (63). 
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Figure 6. (Colour online) The (a) and (b) components of the susceptibility, 
in units of 1/ J - a comparison of the RPA (black solid line) , LSW (blue dotted line) 
and MFA (red dashed line) results below Tn (for d/J = 0.02, AT/ J = 0.42 x IQ-^, 
AJi/J = 0.4 X 10-3). 

Figure 6b shows the z component of the susceptibility x^, and again we obtain the 
good agreement between the RPA and LSW methods at low T, and coincidence of all 
results in the zero temperature limit. On the other hand, in the vicinity of the transition 
temperature, the RPA scheme gives qualitatively different behaviour of with respect 
to the MFA and LSW formalisms. Thus, we can answer one of the motivating questions 
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of this study: Does the extension of the model of reference [3] from 2D to 3D lead to a 
reduction of the strong effects of quantum fluctuations? Indeed, the answer is no, and 
there are strong effects of quantum fluctuations in our 3D Heisenberg model with the 
anisotropics. This statement is correct for a magnitude of the net interplanar coupling 
AJj_ up to ~ 10" V. 

Now let us find the correlation between the ratio of spin- waves modes of the magnon 
excitation spectrum in the long wavelength limit, and the behaviour of the components 
of susceptibility in the zero temperature limit (similar to the correlation that we have 
found in the pure 2D model [3]). Firstly, from the analytical results, equation (63), we 
obtain that the ratio between the components of the susceptibility is given approximately 

by 

- (70) 



T^o 4 J(Ar + 2A Jx) ■ 

Next, by taking into account that the canted angle is ^ ~ (d/-\/2)/(2J) and -C J, we 
can rewrite the expressions (66)- (69) that specify the spin- wave gaps, which we write in 
a scaled form using £i = 2ZrjuJi as 

ul ^ JAJ± + ujI ^ ul ^ J(Ar/2+AJx), ul ^ JAr/2. (71) 

Thus, the ratio between the components of the susceptibility turns out to be 



X 



y 



e2\ 



2 



r->0 \£3/fe-^0 



We also find that gap Ei is always greater than £2, and £3 is always greater than 
£4, when AJ_|_ = J_|_ — > (indeed, that is the only case considered in this paper 
- see the discussion of figure 1 in section 2). As was mentioned above, the two modes 
£1, £2 describe the in- plane modes of the spectrum, while the modes £3, £4 describe the 
out-of-plane spin-wave excitations. Therefore, we find that the observed ratio between 
the X and y components < (in the T — Q limit) [1], in any of the MFA, LSW 
theory, or the RPA, takes place only if the spin-wave gaps have the following hierarchy: 

£1 > £2 > £3 > £4, (73) 

i.e. the in-plane modes (£1,2) are greater than the out-of-planc ones (£3,4). 

The described situation is presented in figure 7 - they show the susceptibility 
for the different values of the interplanar parameter AJ^. The upper curve was 
obtained for the pure 2D case and corresponds to the situation with the observed 
order of the susceptibility components x^ < and the following ordering of the gaps 
£1 = £2 > £3 = £4. By increasing the magnitude of the interplanar parameter AJ^ 
two modes £1 and £3 increase and the hierarchy of the gaps (73) remains unchanged 
(figure 7b). As we can see from the middle curve in figure 7a, the ratio Ix^ decreases 
as the ratio between the gaps £2/^3 decreases. The magnitude of the out-of-plane mode 
£3 becomes equal to the in-plane one £2 when AJ_l/ J fa \{\{d/ J)^ — AF/ J) a; 2.1 x 10~^, 
where the ratio x^/x^ go^s to unity. A further increasing of AJj^ changes the ratio 
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Figure 7. (Colour online) (a) All three components of the susceptibility within the 
RPA analytical scheme for different values of the interplanar anisotropy AJ±: 2D 
result - upper curve, AJ±/J = 0.1 x 10~^ - middle curve, and AJj_/J — 1.0 x 10"^ 
- top curve; and the (b) T — energy gaps vs. the interplanar anisotropy AJ±, for 
d/J = 0.058, AT/J = 0.42 x 10"=^. 



between the modes £2/^3, and according to equation (70) changes the order of the 
susceptibihty components x, z and y at zero temperature (see the lowest curve in 
figure 7a and the corresponding value of the gaps at AJ^ = 0.001 in figure 7b). 

For completeness, in figure 8 we present the susceptibility and behaviour of the 
gaps vs. AJ^ for smaller magnitudes of the in-plane anisotropy. We find that the 
interplanar coupling introduced into the problem leads to a suppression of the 2D 
quantum fiuctuations caused by intraplanar anisotropics. For the large magnitude of 
AJ± = 10^^, i.e. AJ± ~ 2Ar ^ d, the net interplanar coupling AJ± dominates over 
the DM and XY-like pseudo-dipolar anisotropics (see the lower curve in figure 8a). 

Finally, we conclude the presentation of these numerical results by returning to 
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Figure 8. (Colour online) (a) All three components of the susceptibility within the 
RPA analytical scheme for different values of the interplanar anisotropy AJj^: 2D 
result - upper curve, AJ±/J = 0.1 x 10"'^ -middle curve, and AJ±/J — 1.0 x 10"'^ 
- top curve; and the (b) T — energy gaps vs. the interplanar anisotropy AJl, for 
d/J = 0.02, Ar/J = 0.42 X 10"^. 



a discussion of the correlation between the magnitude of the zone-centre spin-wave 
gaps and the behaviour of the components x^'^? X^- Oiig can see that only when 
is greater than 62 at zero interplanar coupling AJ±_ = does the y component of the 
susceptibility at T = become less than components x^'^, since the ratio between gaps 
remains unchanged for any values of AJ± (see equation (72)). 

4. Approximate simple tetragonal model 

For the model parameters of interest our initial Hamiltonian (1) can be approximated 
by a simple tetragonal model Hamiltonian which includes the intraplanar isotropic 
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Heisenberg interaction J, the anisotropic DM term D that alternates in sign from 
bond to bond, the XY-hke pseudo-dipolar anisotropy AF, and an effective interplanar 
interaction J'^/^ . This effective model is defined by 

H = 5^[JS, ■ S, - ^VStS^^ + D,, ■ (S, X S,)] + ^ r/^Sk ■ Sk' . (74) 

{k,k') 

The single-plane effective Hamiltonian was proposed by Peters et al. [16] long ago, and 
its reliability was demonstrated in our previous work [3] . Here the interplanar coupling 
J^-^-^ is added phenomenologically for a simple tetragonal lattice (see figure 9), i and j are 
nearest-neighbour sites in the same CuO plane (indexes ii, ji and 12,32 in figure 9) and 
k, k' are nearest-neighbour sites in adjacent planes (indexes 11,12 and ji, j2 in figure 9). 
Since the interplanar coordination number for interplanar interaction within a simple 
tetragonal model is half that of the corresponding one for the coupling AJ^, we can 
approximate .P/^ = 2AJ±. 




Figure 9. (Colour online) Magnetic structure of the simple tetragonal lattice described 
by the effective model Hamiltonian of equation (74). Use of this Hamiltonian for this 
lattice accurately approximates the magnetic response of the La2Cu04 crystal in the 
LTO phase. 

We performed the calculations of the order parameter, spin-wave excitations and 
susceptibility within the RPA scheme for the effective simple tetragonal model of 
equation (74), and found that the transition temperature, spectrum, and behaviour 
of the order parameter and susceptibility almost identical to that on the initial model of 
equation (1). In figure 10 we show representative data for the susceptibility obtained for 
the initial body-centred orthorhombic model as well as for the effective simple tetragonal 
one with J'^-^ = 2AJ^ - clearly, the agreement between the predictions for these two 
models is excellent, so when studying the magnetic properties of the model (1) on 
3D body-centred orthorhombic lattice one can utilize the effective Hamiltonian of the 
simple tetragonal lattice. Consequently, the magnetism of the La2Cu04 system in the 
LTO phase can be modelled by the Hamiltonian of equation (74). 
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Figure 10. (Colour online) The susceptibility, in units of 1/ J, as a function of T/Tjv. 
We present both results from the RPA calculation on the initial model Hamiltonian 
of equation (1) with AJ±/J — 0.1 x 10~^ and the simple tetragonal model with 
the Hamihonian of equation (74) with J^^^ /J = 0.2 x IQ-^ (for d/J = 0.02 and 
Ar/J = 0.42 X 10^^). The curves from these two models essentially coincide, and no 
differences can seen on this scale. 



5. Conclusions and Discussion 

In this paper we presented a theoretical investigation of the body-centred orthorhombic 
lattice Heisenberg antiferromagnet with in-plane symmetric and anti-symmetric 
anisotropies, and a weak anisotropic AFM interlayer coupling. Our study was focused 
on the role of the different interactions in explaining the magnetic properties of a 
La2Cu04 crystal in the low-temperature orthorhombic (LTO) phase. Due to the 
transition into the orthorhombic phase, the AFM interplanar coupling between nearest- 
neighbour spins in the adjacent Cu02 planes exhibits a small anisotropy. We have 
found that such a small anisotropy plays an important role in magnetic properties 
of the system. In figures 7a, 8a one sees a significant change in the behaviour of 
magnetic susceptibihty as a function of temperature by varying the magnitude of the net 
interplanar coupling AJ^. We also obtained that the (larger) individual superexchange 
interaction between any two nearest-neighbour spins in the adjacent planes does not 
affect the physics of the model (figure 3). 

Our results have shown that in the case of an isotropic interplanar coupling, 
2D quantum fluctuations dominate over the effects of the 3D interaction, and the 
transition to the long-range magnetically ordered state, as well as the behaviour of 
the susceptibility, order parameter and magnon excitation spectrum, are not influenced 
by the interplanar exchange coupling (however, for a 3D model the z component of 
the susceptibility will not diverge, as it does in a 2D model). Thus, in the case of 
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the body-centred lattice model of equation (1) with an isotropic interplanar coupling 
{J'j_ — J±) one can analyze the system using a 2D square lattice model with intraplanar 
anisotropics only. 

We have also obtained that the initial model Hamiltonian (1) can be effectively 
replaced by a simpler one with fewer model parameters, namely by the AFM Heisenberg 
Hamiltonian with DM interaction, XY-hke pseudo-dipolar anisotropy, and an effective 
interplanar interaction J'^^'^ (added phenomenologically for a simple tetragonal lattice). 
Here J^-^ /2 ~ 10~^ describes the small anisotropy of the AFM interplanar couphng in 
the initial system. 

We emphasize an important conclusion that can be drawn from our results. We 
have found that in-plane anisotropy introduced into the problem by symmetric XY-like 
pseudo-dipolar and antisymmetric DM interactions largely determines the behaviour 
of the magnetic susceptibility, transition temperature into the long-range order state, 
and the spin-wave gaps in the case of a 3D model (within the wide range of model 
parameters of interest). Further, even when one studies a 3D model, the effect of 
quantum fluctuations is very strong in all temperature regions below the transition 
temperature, and cannot be ignored. Similar to the results of our previous paper [3], we 
have also obtained the large short-range correlations in the broad temperature region 
above the Neel temperature. 

Now we comment on the comparison of our results to the experimentally observed 
anisotropics of the susceptibility [1] and spin- wave gaps [15] that motivated our work. 
We can state that all anisotropic interactions involved in the model, i.e. DM, XY-like 
pseudo-dipolar, and interplanar ones, are responsible for the unusual anisotropy in the 
magnetic susceptibility, and the appearance of gaps in spin-wave excitation spectrum. 
More concretely, by comparing to a purely 2D model, the inclusion of interplanar 
anisotropy leads to the splitting of either of the in- and out-of plane zone-centre spin- 
wave modes. While the neutron-scattering experiments find only two gaps, one in-plane 
mode Ei ~ 2.3 mcV and one out-of plane mode ~ 5 meV, we can infer the following 
possible situation that is predicted from our results: the in-plane mode ei (which is 
always larger than 82) has a gap with a magnitude of about 10 meV. Indeed, such an in- 
plane gap can be seen from the result for the spin- wave spectra in the neutron-scattering 
experiments [15]; other observed gaps corresponds to the out-of plane mode £3 ~ 5 meV 
and in-plane mode £2 ~ 2.3 meV. The magnitude of the gap of the remaining out- 
of plane mode, £4, is relatively small and apparently has not be seen by experiment. 
Therefore the following hierarchy £1 > £3 > £2 > £4 agrees with the experiment. In this 
paper we established the correlation between the ratio of the in and out-of-plane gaps of 
the excitation spectrum and the behaviour of x^'^ and components of susceptibility in 
the zero temperature limit. However, the proposed hierarchy of the spin- wave gaps takes 
place only if the ratio between the x and y components is opposite to that observed in 
experiment (x^ < x^)- This necessarily leads to the question, would other interactions, 
e.g. ring exchange and/or the interaction between the next nearest neighbour sites [6], 
lead to an accurate explanation of the susceptibility data within the RPA scheme? 
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In order to answer this question we have performed calculation for the square 
lattice AFM Heisenberg model with the DM and XY-like pseudo-dipolar anisotropies 
by additionally taking into account the ring exchange and the interactions between the 
next nearest neighbour sites (for the energy scales of these additional interactions see 
reference [6]). Our results of the RPA calculations have established that ring exchange 
together with the second and third nearest neighbour in-plane exchanges do not change 
the results presented in our earlier paper [3] regarding the correlation between the ratio 
of the in and out-of-plane spin-wave modes gaps and the behaviour of the 
components of susceptibihty in the zero-T limit, viz. el/ef ^ X^/x^- So, physics beyond 
what has been presented in our previous and this manuscript is important, but that does 
not imply that a more complicated Hamiltonian with more interactions is necessarily 
required. 

A potential resolution of this dilemma can be found in studies based on the 
quantum non-linear sigma model [14] . However, within a theory that accounts for short- 
wavelength behaviour we will show in a future publication how the "next" approximation 
beyond that used in our previous [3] and present papers fits the experimental data. This 
allows for the important next problem of the couphng of the anisotropic AFM state to 
either localized or mobile holes to be examined. 
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The transformation between the initial representation and the characteristic 
representation (CR) in which the quantization axis is in the direction of the classical 
moment (see figure 1) reads 
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In our calculations we formally divided the body-centred lattice into the four 
sublattices which differ in orientation of their spins within the ground state. As was 
explained in section 2, we follow the notation where ii-sites belong to sublattice 1, 
ji-sites to sublattice 2, i2-sites to sublattice 3, and j2-sites to sublattice 4, respectively. 

By using the CR as the basic representation we obtain the components of 
susceptibility that determine the response of the expectation value of the spins in 
one sublattice with respect to the external magnetic field applied to another one. For 
instance, the component 



^ iV/4 7V/4 
A 



EE [{TralirKmdr , (A.5) 



ii=l ii=li2=l' 

determines the response of the expectation value of the spin a^, 4/A ^^^''(crf^), in 
sublattice '1' to the field applied to sublattice '3' in the cr^-direction of the corresponding 
coordinate system within the CR [3] . 

The transformation between the susceptibility components in the initial and CR 
coordinates reads 

^.x f^.o-'^cr^ _i_ ^.a-'^a'^ i ^.a^a^ i ^.a'^a^ -v/f'^f^ -v/f"^"'^ -v/f^f^ -v/f"^^^ 

A — 2 lAll "T A12 "T A13 "T A14 All Al2 A13 A14 



"T All "T A12 "t" A13 "T A14 All Al2 A13 A14 J ' 

A — lAll "T A12 ~ A13 ~ A14 All "r A12 ~ A13 ~ A14 



(A.6) 



"T All "T A12 ~ A13 ~ A14 All "r A12 ~ A13 ~ A14 J 

+ cos'{9){xr^''' - Xi2^^ - Xn"^ + xf/l, 



(A.7) 
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Appendix B. System of equations for the Green's functions in the 
momentum- frequency representation 

Let us present the most general form of the system of equations for the Green's function's 
in equation (25). We introduce notation g for the zero-order gin = Gin{k,uJm) and/or 
first-order gin = G[1^ {k,ujm) Green's functions depending on the coefficients in the 
equations (see below). Then, the system reads 

^^9i2 = Ui2 + -^9i2 + Akg22 - -Bfe5'22 + afe5'42 - iOfe5'42 - l"fe5'32 > 

Jmfa 



„ ^ 5-22 = t^22 + —^922 + Akgi2 " 3^912 + akg32 - ^^932 - ^dk942 > 

i^m — j-T— Jmfa _ . _ . 7-,^= _ . , . , 

^^fl'12 = - f^l2 ^9l2 - ^k922 + ^fefl'22 - afc5'42 " lOfcfi'42 " iafc5'32 , 

fl'22 = - ^^2^2 - -^922 - '^k9i2 + B*k9i2 - ak932 - i^fefl'32 - ic^fe5'42 , (B.l) 



icUm _ Jmfa — 

i^m TT- . Jmfa 



22^332^ U32 + —^932 +Ak942 + Blg42 + (^k922 -ibk922 -^dk9l2 > 
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i'-'-'m _ Jmfa _ 4 _ t-, _ ., . , 

^^942 = - (^42 2~^^^ ~ '^k932 " -Bfefl'32 " afc5'l2 " Wk9l2 " iafe5'22 ■ 

In the case of the zero-order system, gin = Gin{k,u!m), we have 

U22 = Ui2 = Ur2 = C/2"2 = t^32 = C/42 = U^2 = t^4"2 = 0. (B.2) 

In case of the first-order system gin = G\l^{k, Um) we have 
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where the new quantities Vin and were introduced. 

Appendix C. Zero-order Green's functions 

For ease of presentation, we define some new quantities: 

xi,k = -2ak[ak{Ak-Jmfa/'2) - (6fc-4)^>5fc], yi,k = -(^fc- Jm/a/2), 

X2,k^ '^ak[ak{Ak+Jmfa/'2) - {bk+dk)'^Bk], y2,fe = (^fc+J^/a/2), 

X3,fe = 23?Sfc[-afc + (6fc-(ifc)^] + 2i{bk-dk)[ ak{Ak-Jmfal'i) - '^Bk{bk-dk)], 
X4,k = 23fJ5fc[ al - {bk+dkf] + 2i{bk+dk)[-ak{Ak+Jmfa/2) + '^Bkibk+dk)], 

y3,k = -Bl, 
y4,k = Bi, 

X5,k = -2{Ak-Jmfa/'2)[ak{Ak-Jmfa/'2) " (6fe-dfc)3^Sfe] + 2afc3?^Sfe, 

+ 2mBk[(Ak-Jmfa/'2)(bk-dk) - akQBk], 
xe^k = 2(Ak+Jmfa/'2)[ak(Ak+Jmfa/'2) - (bk+dk)QBk\ - 2ak^^Bk, 

- 2mBk[{Ak+Jmfa/'2){bk+dk) - ak^Bk], 
y5,k = -Qfe, ^5,fe = -2Q^-Bfe(&fe-dfe) + 2ak{Ak-Jmfa/'^) - 2i3?-Bfe(&fe-dfc), 
?/6,fe = Ofc, ze^k = -2QBk{bk+dk) + 2ak{Ak+Jmfa/'2) - 2i?R:Bk{bk+dk), 

x-r^k = 2i[akiAk-Jmfa/'^)^Bk - ibk-dk)\Bk\'^], 
X8,k = -2i[ak{Ak+Jmfa/'2)^Bk - {bk+dk)\Bk\'^], 
y7,k = i{bk-dk), zr^k = -2ak^Bk, 
ys,k = -i(&fe+o?fe), 2;8,fe = -2ak^Bk. 

Then, the solution of the system of equations for the zero-order Green's functions can 
be written as 
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(C.3) 



(C.4) 
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I^S.fc ^ i^m — £3,fc il^m + ^3,fc ^ 
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(C.8) 

Appendix D. First-order Green's functions 

Prom the whole set of the first-order Green's functions we need only the diagonal G^'^ 
ones. The solution for such Green's functions can be presented via the coefficients 
(defined in equation (B.3)) and zero-order Green's functions in the following form: 

Gff = -Z { \Gu\'V22 + 1^221 Vi2 + |G'r2| V2i + IG2-2I (D.l) 

+ 1^32^142 + 1^421^1^32 + 1^321^142 + l^42r^2 } 

G22 — ~Z I |Gl2|^1^12 + |G^22|^1^2 + I^12r^2 l^22l^^2 (■'-^■2) 

+ |G'32|^V32 -|- |G'42|^V42 -|- 1^321^132 + l^42l^Kl2 } 

= -Z { IG12IV42 + IG22I V32 + |Gr2l + IG2-2I V32 (D.3) 

+ |G'32|^V22 + |G42|^Vi2 + |G'32|^V22 + 1^421^142 } 
= -Z { 1^121^32 + 1^221^42 + \G^2\%'2 + |G'2-2l'V;2 (D.4) 
+ 1^321^^12 + |G42|^V22 + |G^32|^^12 + 1^421^^22 } 
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Appendix E. Transverse components of the susceptibility in the 
characteristic representation 



We can find the transverse components of susceptibility from the following relations 

(E.l) 



for any sublattice /,n = 1,2,3,4. The components " i X" can be found 
immediately from the solution of the zero-order Green's functions in Appendix C and 
the definition in equation (50) with a = +, — : 
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where all coefficients xi, yi with / = 1, . . . 8, which are taken from the Appendix C, and 
all frequencies uj, and /5i^2 (see equations (41,43)) are taken in the long wavelength limit 
fe — > 0. Note that in above formulae we have also used the following relations between 
the different components of the transverse components of susceptibility in the CR 

Xll = X22 = X33 = X44, Xl2 = X21 = X34 = X43, (E-ll) 

Xl3 = XSI = X42 = X24, Xl4 = X41 = X32 = X23- 



Magnetic susceptibility of La2CuOi 



34 



References 

Lavrov A N, Ando Y, Komiya S and Tsukada I 2001 Phys. Rev. Lett., 87 0170071 

Tranquada J M 2005 cond-mat/0508272 

Tabunshchyk K V and Gooding R J 2005 Phys. Rev. B, 71 214418 
Dzialoshinski I 1958 J. Phys. Chem. Solids, 4 241 
Moriya T 1960 Phys. Rev., 120 91 

Katanin A A and Kampf A P 2002 Phys. Rev. B 66, 100403(R) 

Johnston David C 1997 Handbook of Magnetic Materials 10 (Elsevier, New York, 1997) 
Coffey D, Rice T M and Zhang F C 1991 Phys. Rev. B, 44 10112 
Shekhtman L, Ebtin-Wohlman O and Aharony A 1993 Phys. Rev. Lett, 71 468 
Koshibae W, Ohta Y and Mackawa S 1994 Phys. Rev. B, 50 3767 

Xue W, Grest G S, Cohen M H, Sinha S K and Soukouhs C 1988 Phys. Rev. B, 38 6868 
Lee K H, Liu S H 1967 Phys. Rev., 159 390 
TyabHkov S V 1959 Ukrain. Mat. Zh, 11 287 

Silva Neto M B, Bcnfatto L, Juricic V and Morais Smith C 2005 cond-mat/0502588 
Keimer B, Birgeneau R J, Cassanho A, Endoh Y, Greven M, Kastner M A and Shirane G 1993 Z. 
Phys. B, 91 373 

[16] Peters C J, Birgeneau R J, Kastner M A, Yoshizawa H, Endoh Y, Tranquada J, Shirane G, Hidaka 
Y, Oda M, Suzuki M and Murakami T 1988 Phys. Rev. B, 37 9761 



